Universal correction for the Becke-Johnson exchange potential 
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Abstract 

The Becke-Johnson exchange potential [J. Chem. Phys. 124, 221101 (2006)] has been success- 
fuhy used in electronic structure calculations within density-functional theory. However, in its orig- 
inal form the potential may dramatically fail in systems with non-Coulombic external potentials, 
or in the presence of external magnetic or electric fields. Here, we provide a system-independent 
correction to the Becke-Johnson approximation by (i) enforcing its gauge-invariance and (ii) mak- 
ing it exact for any single-electron system. The resulting approximation is then better designed to 
deal with current-carrying states, and recovers the correct asymptotic behavior for systems with 
any number of electrons. Tests of the resulting corrected exchange potential show very good results 
for a Hydrogen chain in an electric field and for a four-electron harmonium in a magnetic field. 
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Electronic structure calculations are routinely carried out by using density-functional 
theory |2[ (DFT) and its variants. The accuracy of a DFT calculation depends on the 
approximation used for the exchange-correlation energy functional. Substantial efforts have 
been made in deriving accurate approximations over the past few decades j^. 

Within spin-DFT (SDFT) the optimized-effective-potential (OEP) method jj-lsl provides 
a rigorous access to the exact exchange (EXX) energy, 



(1) 



a-=t,ij,k=l ■ 

and to the Kohn-Sham (KS) exchange potential fa,.(T(i") = SE^/ Pair). Hartree atomic units 
(a.u.) are used throughout. Above v^jo-(r) are the spin-dependent Kohn-Sham (KS) orbitals, 
with energies Eja, and 

is the ground-state density. The OEP method leads to an integral equation which can be 
solved iteratively together with the standard KS equations. The main origin of practical 



complexity in the OEP method are the orbital shifts containing uno ccup ied KS orbitals 



Despite great pro gres s in solving the equations for various systems [sHlGj , and in algorithmic 
developments 17|, ll8| , efficient and accurate approximations reducing the numerical burden 
of the full OEP scheme are still needed. 

A simple approximation for the exchange potential has been proposed by Becke and 
Johnson (BJ) 
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is the Slater potential, 

is (twice) the spin-dependent kinetic-energy density, and Cav = [5/{12n^)]'^^^ . Note that the 
Slater potential in Eq. @, which is now the only numerical bottleneck, could be alterna- 



tively approximated by the semi-local Becke- Roussel approach 



20|. Interestingly, the exact 



exchange potential could be written as in Eq. (|3]), but with Af^^(r) replaced by Aw^f^i 
which can be decomposed into the so-called Krieger-Li-Iafrate (KLI) approximation 21| . 

n 

plus another correction given in terms of the orbital shifts [7]. As it is well known, in many 
cases the KLI approximation is in good agreement with the full OEP. 

Despite the semi- locality of Af^^(r) in Eq. ([2]), the BJ potential is able to correctly 
yield the step structure in the exchange potential of several atoms [l9|. Moreover, it has 
recently been shown that the BJ potential correctly reproduces the derivative discontinuity 



"ter its introduction, the BJ 
28l |. Impressively, the band 



for fractional particle numbers [2^ . During the first few years a: 
approximation has already been applied to various systems 22- 
gap of a large variety of extended systems is extremely well reproduced [26|, |28|] . However 
as we will demonstrate below, the BJ potential may dramatically fail in the presence of an 
electric or magnetic field, or a non-Coulombic external potential. 

The limitation of the BJ potential originates from two facts: (i) it is not gauge-invariant 
and (ii) it is not exact for all possible one-electron systems. These two problems may be 
fixed in similar fashion as demonstrated in our recent derivation of a BJ-type approximation 
for two-dimensional systems 29 1. 

Before proceeding further, we would like to comment on the gauge- invariance requirement. 
For systems acted upon an external vector potential, the exchange potential we propose in 
this work should be identified as an approximation derived for the exact exchange-potential 



obtained within current-spin-density functional theory 30|, |3l| (CSDFT) by taking the func- 
tional derivative of the exchange-correlation energy functional (written in terms of the spin- 
particle and vorticity density) at constant vorticity. On the other hand, it is clear that since 
Ex in SDFT depends only on Po-(r), it must be a gauge- invariant quantity by definition. As 
a direct consequence, the corresponding fa;o-(r) is gauge-invariant as well. 
Therefore, we propose the following correction 
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where 



j^-(^) = ^ E (V^..(r)) - (V^*,(r)) v,,.(r)] 



(8) 



is the spin-dependent paramagnetic current density. The above potential in Eq. (|6]) has a 
set of desirable properties listed below. 

• In contrast with Ta-(r), as it appears in the BJ expression, the combination T„{r) — 
jp^(r)/po-(r) is clearly gauge- invariant. |32M34| As a result, also the corresponding 
potential is gauge-invariant. 



In contrast with Aw^^, Af^^ is zero for all one-particle systems. This is easy to see 
by considering an arbitrary one-particle system with Po-(r) = |y9o-(r)p and v^cr(r) = 
^y^^^^^^gje(r)^ so that D„{y) = follows immediately from Eq. ([7]). Alternatively, this 
may been seen by using the definition of To-(r) and jpo-(r) [Eqs. ([5]) and (|8]), respectively] 
in terms of V'o-(r), and by re-expressing Po-(r) = |(y9o-(r)p in the second term of Af£^. 

The asymptotic limit is correct for any A^-electron finite system: Af2^(r — )■ oo) — ?■ 
and then Vxaij oo) — > v^{v — )■ oo) — )■ — 1/r. In that limit all the terms in 
are dominated by the highest occupied orbital [35], and thus the system effectively 
behaves like a one-particle system (see the preceding point). Below we discuss the 
asymptotic limit in detail for two particular systems. 

Equation (E]) is consistent with the limit of the homogeneous 3D electron gas (3DEG): 
At;O^P = ^vf^ = Avg = [3p./(4vr)]^/l 



Calculation of Af£^ instead of Av^^ does not bring any extra computational burden. 



Also, it is reassuring to note that, the exchange potential in Eq. scales linearly as 
the exact one (see Appendix) [36 1. 



Finally we point out that the key object in the corrected exchange potential, Dcr{r), is 
familiar from various concepts in the literature. First, it is an important part of the electron 

37H39|. and second, it enters in the expression of the local curvature o: 



localization function 



the exchange hole j40|. In the latter case, it is a part of the current-generalized forms |4l|-|43| 
of the Becke-Roussel and Becke models for the exchange [2^ and correlation 44 1, respec- 
tively. 
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FIG. 1: (Color online) Difference in the exchange potentials for a four atom hydrogen chain with 
and without the linear electric field. Black dots denote the positions of Hydrogen atoms. 



Next we test our exchange potential against the KLI [21], BJ, and local-density approx- 
imation (LDA) for two different systems. We perform the self-consistent KLI calculations 
applying the octopus j45| DFT code on a real-space grid. The resulting KS orbitals are 
then used as inputs in the approximations for the exchange potentials. 

First we consider a H4 chain in an external linear field with the same system parameters as 
in the work by Armiento, Kiimmel, and Korzdorfer 22| (AKK). The system consists of two 
H2 "molecules" with an interpair distance of 2 a.u. separated by 3 a.u. The strength of the 
electric field, applied along the x direction, is F = 0.005 a.u. (hartree/bohr). Figured] shows 
the difference in the exchange potential with and without the electric field, i.e., v^aiF = 
0.005) — Vxa{.F = 0). Overall, we find excellent agreement between KLI (solid line) and the 
present approximation (dashed line) in Eq. ([6]). 



Figure [T] can be directly compared to Figs. 3 and 5 in Ref. 22|. We find the same 
divergence in the BJ potential in the asymptotic regime, as well as the large deviation of the 
LDA from the KLI result. More importantly, we find that our corrected formula is closer to 
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KLI as the potential proposed by AKK, i.e., 




X F 



-2e^^^ - ^=== I , (9) 



with en^o- corresponding to the energy of the highest occupied KS orbital |46[. In particular, 
the AKK potential difference is close to zero on the left (x < —10), whereas the KLI and the 
present approximation yield a finite value in that regime. Close to the atoms we find some 
overestimation in the maxima (and minima) of our potential, but, on the other hand, our 
approximation is free from sharp kinks present in the BJ and AKK potential differences. 

Next we compare in detail the asymptotic limit of the above example given by the AKK 
potential [Eq. ([9])] and our potential [Eq. (E])], respectively. For atomic systems in the 
absence of external fields the asymptotic behavior of the KS orbitals is given by Q 



(^,.(r) ^ $,.(r)/,.(^]) , (10) 

with $jo-(^) being the asymptotic radial wavefunction, and fja{^) its corresponding angular 
component. The asymptotic form of $jo-('") is jj] 

(11) 



with (3 jo- = 1^—2 Ejo-. Substituting Eqs. (fTOj) and (fTT]) into Eq. ([6]), it is easy to find that 
the leading correction in the asymptotic limit is given by 



.1/2 



(12) 



where gN^a{^) = - ["^ fN„a{^)]['^ fN^oi^)]/\fN^a{^)\'^ is a purely angular term with con- 
tributions coming from the second and third terms in defined in Eq. ([7]). Making the 
subsequent asymptotic expansion of r^ir) / po{r) '"^°°> —2eN^cr— gN^ai^), it is apparent that 
the term inside the square-root in Eq. f|T2|) vanishes identically in the asymptotic regime. 
We emphasize that, Eq. (|T2l) is similar, but not identical, to the AKK potential in Eq. (j9]), 
with F = 0. The main difference is that while the satisfaction of several exact constraints 
(as explained above) enforces us to have all the contributions of inside the square-root, 
the correction in the AKK potential that enforces the vanishing of the exchange potential 
in the asymptotic limit {y/—2eN^) is outside the square- root. 



Similar considerations apply when the external electric field F is present. Following AKK, 
in this case the asymptotic behavior of the KS atomic orbitals along the direction of the 
applied electric field is given by jl^] 

with 1] = {2Fyf^{x — Eja/ F). Proceeding with the evaluation of At>^^(x — )■ oo), we find 



2{xF - En^^) 



1/2 



(14) 



which again is similar, but not identical, to Eq. (Q. Making them identical would require 
an expansion of the argument inside the square-root, but this is unjustified, since in the 
asymptotic regime both contributions are equally important. That is, t„{x) / pa{x) ^~^°°> 
2{xF — EN^a), leading to a cancellation of both terms inside the square-root. 

As a conclusion of the analysis in two previous paragraphs, let us emphasize that both our 
f3;o-(r) and f^J^^(r) reproduce the correct asymptotic limit of the exact exchange potential, 
but in different ways. On the other side, v^^^{r) is system-dependent, being only valid for 
atomic systems in presence of a bias, while our exchange potential is system- independent, 
being valid for any 3D system, in the presence of any electric and/or magnetic fields (see 
below). Also it is better suited for dealing with currents. Moreover, our potential does 
not require explicit knowledge of external fields (as it should be for any standard density 
functional) and/or KS eigenenergies - only (occupied) KS orbitals are needed. 

In Fig. |2] we show the exchange potentials for a fully spin-polarized four-electron "har- 
monium", i.e., a 3D harmonic oscillator with a radial external potential fext('") = 
where a; = 1 a.u. This type of potential could be used as a realistic model for quantum 
dots, i.e., electrons confined in atomic clusters or semiconductor heterostructures. We have 
also set an external, uniform magnetic field to -B = 300 a.u., so that the occupied KS states 
have angular momenta / = 0, — 1, —2, —3 and hence there are orbital currents in the system. 
The BJ potential (thick dashed line) shows erroneous divergent behavior. Similar divergence 



appearing in two-dimensional harmonic oscillator has been analyzed in detail in Ref. 29 1. 
It was shown that the linear increase in the BJ potential at large r follows directly from 
the asymptotic limit of the single-particle wave functions, which, in the case of a parabolic 
confining potential, decays as exp(— r^) in contrast with the atomic wave function that de- 
cays as exp(— r). The situation is the same in the 3D case considered here. The LDA result 
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FIG. 2: (Color online) Exchange potentials for a four-electron harmonium in external magnetic 
field B = 300 a.u. 

(dash-dotted line in Fig. [2]), on the other hand, largely underestimates the exchange poten- 
tial throughout the system. Similar tendency is shown by the BJ potential modified by the 
gradient term [second term in Eq. ([71)] but without the current term [third term in Eq. (^^], 
i.e., without enforcing the gauge-invariance (thin dashed line marked by v^^ + At>°^). The 
closest resemblance of the KLI potential in Fig. [2] is clearly given by the present approx- 
imation in Eq. ([6]). Also it can observed that although the asymptotic limit is very well 
reproduced, close to the core of this system we still find some deviation, both for the BJ and 
for our corrected exchange potential. This gives evidence that further improvements may 
be suggested in future works. 

Finally, we verify that the atomic step structure at electronic shells - one of the moti- 



vations behind the original Becke- Johnson approximation [19| - is reproduced by the cor- 
rected potential. In Fig. |3] we show the exchange potential for a closed-shell neon atom in 
its ground state with no external fields present. The corrected potential (dashed line) has 
the step structure at r ~ 0.3 a.u. in agreement with the BJ potential (dotted line) and 
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FIG. 3: (Color online) Exchange potential for a closed-shell neon atom in the ground state. 



with the KLI potential (solid line). The OEP solution shown in Ref. 19| shows a slightly 



sharper shoulder than the KLI one, but generally they are very similar, which is in ac- 
cordance with previous studies in the literature (see, e.g., Ref. |7|). In the exchange-LDA 
potential (dash-dotted line) the step (or shoulder) structure is missing. Note that for con- 
sistency with the previous results we have not iniposed a shift to the BJ potential, which 



is a prerequisite having different definitions (l9|, |22|. Secondly, we point out that Fig. [3] 



results from self-consistent calculations for all potentials, respectively. The difference from 
the non-self-consistent results, obtained by using the KLI orbitals as the input, was found 
to be negligible. 

We point out that for atoms at small r our potential decreases more strongly than the 
KLI (or BJ) potential (see Fig. [3]). In fact, at r = the correction term Af^ is significantly 
smaller than Av^^^. This difference is not present at r = in a harmonic confinement 
discussed above, since in that case it is easy to show that the correction terms in the BJ 
and our potential have exactly the same value due to the Gaussian form of the orbitals; 
this feature can also be observed directly from Fig. [2l Nevertheless, the deviation found 
in the atomic case close to the nuclei suggests that an additional effort beyond the present 
contribution may be required. 
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Concluding, we have seen that Becke- Johnson potential may dramatically fail when con- 
sidering systems in external electric and magnetic fields. We have proposed a universal 
correction which is gauge-invariant for complex Kohn-Sham orbitals and exact for any one- 
particle system. The improved approximation is suited for dealing with with current-carrying 
states, and it also recovers the correct asymptotic behavior of the exact exchange potential 
for any many-electron system. We have demonstrated the very good performance of the 
resulting exchange potential by considering a hydrogen chain in an external electric field as 
well as a four-electron harmonic oscillator in a magnetic field. 
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Appendix 



We may write the term beyond the Slater contribution to the exchange potential as 

'DJr 



Cav 



(A.l) 



and determine a under the constraint of exact linear scaling 36|. Under uniform scal- 
ing of the coordinates, r — )■ Ar, the norm-preserving many-body wavefunction is given 
by \E'A(ri, ...Fat) = A^^/^\['(Ari, Afat) (with < A < oo). As a consequence, the 3D 
density scales with A as follows: Pcr(r) — )■ \'^pa{\r). This leads to the result that the 
KS orbitals in 3D are seen to scale as v'jo-(r) — t- A^/^ (y9jo-(Ar). Thus, ra-(r) — t- A^ro-(Ar), 
Vpcr(r) — > A'' VArP(T(Ar), and jpo-(r) — ?• A^jpo-(Ar). Substituting these relations into Eq. ( lA.ip 
yields Av^J'i^r) = A^"At>£^(Ar), which fulfills the linear scaling constraint only if a = 1/2 in 
agreement with the expression in Eq. (|6]). 
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